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We present Monte Carlo estimates for site and bond percolation thresholds in simple hypercubic 
lattices with 4 to 13 dimensions. For d < 6 they are preliminary, for d > 6 they are between 20 to 
10 4 times more precise than the best previous estimates. This was achieved by three ingredients: 
(i) simple and fast hashing which allowed us to simulate clusters of millions of sites on computers 
with less than 500 MB memory; (ii) a histogram method which allowed us to obtain information for 
several p values from a single simulation; and (iii) a variance reduction technique which is especially 
efficient at high dimensions where it reduces error bars by a factor up to ~ 30 and more. Based on 
these data we propose a new scaling law for finite cluster size corrections. 

PACS numbers: 05.10-a, 64.60.A 



In spite of decades of intensive studies jjj , percolation 
remains an active subject of research. While there has 
been enormous progress in understanding percolation in 
2 dimensions [^J, mainly because of conformal invariance, 
progress in high dimensions has been much slower. It has 
been known since long time that d = 6 is the upper criti- 
cal dimension jlj], and expansions of p c in l/(2d — 1) have 
been given already more than 20 years ago. But up to 
now there exists no detailed numerical study of logarith- 
mic corrections in d = 6, finite size corrections are not 
yet understood for d > 6, and even numerical estimates 
of p c in > 6 dimensions are very poor. One reason for 
this is obviously that straightforward simulations of large 
lattices in > 6 dimensions require huge amounts of fast 
memory. This lack of stimulus by numerical verification 
certainly was part of the reason for the slow analytical 
progress. 

It is the purpose of the present note to improve on this 
situation by presenting precise numerical estimates of p c 
(and of finite cluster size corrections) for site and bond 
percolation on simple hypercubic lattices with d = 6 to 
d=13. 

Our main results are summarized in Table 1, where 
we also include preliminary results for d — 4 and d = 
5. We also give the best previous estimates for p c and 
expansions in 1 / (2d — 1). We shall discuss them later 
in more detail, but here we just point out that our new 
estimates are vastly better than all previous ones. They 
were possible, with rather modest effort (we used only 
fast PCs and Alpha work stations, with altogether ca. 
10 3 CPU hours), due to several important ingredients: 

1) We used as basic routine a standard breadth- first 
version of Leath's algorithm which simulates single clus- 
ters. We do not use the popular Hoshen-Kopelman 
method since that would require prohibitively large mem- 
ory if we want to simulate large clusters. In Leath's 
method, one writes the coordinates of each cluster site 
(which consist of a single integer - see item 2 below) into 
a first-in- first-out queue Q, where each new entry repre- 
sents a newly wetted neighbor of the oldest entry in the 
queue. 

2) We used a simple but very efficient form of hashing 



H for storing the information whether a site has already 
been wetted or not. On Compaq Alpha work stations 
with 64 bit long integers, we labelled lattice sites by a 
single long integer. Using as lattice size L an odd number 
slightly smaller than 2 64 / d , we label the neighbors of site 
i as i± 1, i± L, . . .i± L d_1 . If we want to simulate bond 
percolation clusters with roughly N sites, we find first the 
power of 2 nearest to N, 2 k ss N, and use it to obtain for 
each site i its key = (i mod 2 k ) (notice that this is 
done most efficiently by bitwise AND). Assume now that 
site i with key rrii is the n-th site wetted. Then an entry 
is written into the m-i-th element of an array of pointers 
S of size 2 . This element points to the n-th element 
of a structure (L,Q) where Q is the above queue and L 
is a linked list. In Q, the coordinate i is stored. The 
element of L remains empty, if the key rrii had not been 
encountered before. Otherwise, if some other site j with 



the same key m. 



had been wetted in an earlier step 



nf < n, the old element of S (which had pointed to n') 
is written in the n-th element of L. In this way we can 
deal with virtual lattices of 2 64 sites, using 2 k + 2N maK 
storage places, where -/V max is an upper bound on the 
size of clusters to be simulated. The algorithm is slightly 
different for site percolation where a tested site has to be 
excluded from further growth even if it is not wetted, in 
contrast to bond percolation. It also has to be modified 
on machines with only 32 bit long integers where a pair 
of numbers replaces i and a pair of co-prime odd numbers 
L\ and L2, both slightly smaller than 2 32 / d , replaces L. 

This is not as storage efficient as the recent algorithm 
of But it works with usual (pseudo-)random num- 
ber generators (we used the 4-tap generator with period 
29689 _ 2 f gj^ ^ile the algorithm of j| needs a ran- 
dom function generator. The most easily available ran- 
dom function generator today is the Data Encryption 
Standard |6| which is rather slow when implemented in 
software and of unproven quality for this application (it 
was developed for entirely different purposes, and lacks 
any published theoretical justification). 

3) In order to estimate cluster statistics for several 
values of p from a single run at nominal value pa, we 
use a trick similar to the histogram methods used by 
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.1601314(13) 
.118172(1) 
.0942019(6) 
.0786752(3) 
.06770839(7) 
.05949601(5) 
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.04018762(1) 
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.11664888 
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.04018504 
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.1968861(14) 
.1407966(15) 
.109017(2) 
.0889511(9) 
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.04271508(8) 
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.19304456 
.13793629 
.10754047 
.08823220 
.07485431 
.06502556 
.05749265 
.05153203 
.04669616 
.04269312 
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.19880605 
.14004471 
.10848530 
.08871655 
.07512757 
.06519119 
.05759880 
.05160316 
.04674559 
.04272853 



TABLE I: Estimates of p c for bond and site percolation in d = 4 to d — 13. Numbers in round brackets are single standard 
deviations, square brackets refer to the citations at the end of the paper. For d > 9 the best previous estimates backed by 
theory were given by the (presumably asymptotic) expansions (Q) and (^), while Eq.(^) was a heuristic guess. The estimates 
for d = 4 and d = 5 are preliminary, since we do not yet understand the important corrections to scaling in these cases (all 
error bars in this paper include plausible worst case estimates of systematic errors). 



Dickman jll| for the contact process (see also JT^]). If a 
cluster with n wetted sites and b non-wetted boundary 
sites was generated with nominal value po, it contributes 
to the ensemble with p replaced by p with weight 

w = ( P /p n(i~p)/(i~p a )) b . (i) 

Instead of collecting histograms for cluster numbers with 
fixed n and b (which would have led to excessively large 
arrays) we calculated on the fly three distributions: One 
for the nominal po (which was chosen close to p c as esti- 
mated from short test runs and from Eqs.(Q) resp. (||)) 
and two for neighbouring values p± = po ± dp, using 
Eq.(l) for the latter. Observables at p- values in between 
(including p c ) were obtained by geometric (i.e. linear in 
logarithm) interpolation. Having three values of p in- 
stead of just two allowed us to check that the error due 
to the interpolation was negligible. 

4) Our main observable will be the number M(t) of 
wetted sites with "chemical distance'' t from the seed of 
the cluster (i.e. the number of sites infected at time t, 
if cluster growth is interpreted as spreading of an epi- 
demic). For d > 6 we expect its average (M(t)) to be- 
come a constant at the critical point, since the process 
is basically a branching process with small corrections. 
But instead of using (M(t)) itself, we obtain a less noisy 
signal by the following trick which would give the exact 
ensemble average of M(t) if the cluster growth indeed 
were a branching process Oj. 

Let us assume we have a (still growing) cluster C with 
M(t) sites wetted at step t, and denote by M + (t) the 
number of free neighbors, i.e. the number of sites which 
can be wetted at step t+1. The actual number wetted will 
fluctuate, but the expected average number, conditioned 
on C and thus also on M(t), is exactly given by 

E[M{t + l)\C] = P M+(t). (2) 



Thus the expected geometric increase of the number of 
wetted sites, still conditioned on C, is 

E[M(t + l)/M{t)\C] = P M+(t)/M(t) (3) 

and its weighted sample average over all clusters is 

(f) _ E C M(m[M(t+l)/M(t)\C] _ p(M+(t)) 
{) ~ J2c M (t) (M(t)) ■ [ ! 

Our estimate for the true ensemble average of M(t) is 
then finally 

t— l 

M(t) = J] (5) 

t'=0 

Since we measured also the direct estimate (M{t)) and 
the (co-) variances of both estimates, we can also com- 
pute the variance of any linear combination of both. For 
d = 4 and d = 5, where both variances are comparable 
and the covariance is negative, a substantial achievement 
is obtained by taking as the final estimate the linear com- 
bination with the smallest variance. 

An expression similar to Eq.(||) can be obtained also 
for the rms. radius, if we replace the ratios in Eqs.(||^) 
by differences and the product in Eq. (|5j) by a sum. This 
would also be exact and non-fluctuating if the cluster 
growth were a branching process with translation invari- 
ance. 

The variance reduction due to Eq.(^) is largest for 
small t. Yet, for bond percolation in d = 11, it gave 
even for the largest t (=200) a factor « 1/1000 over us- 
ing just (M(t)}. For d = 6 and t = 2000, the reduction 
was still by a factor « 140. Indeed, there were substantial 
improvements even for d = 4 and 5, while the improve- 
ment in d = 3 was marginal. For site percolation the 
improvements were similar but somewhat less dramatic. 
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FIG. 1: Plot of M(t) versus lnt for bond percolation in d = 6. 
Statistical errors are smaller than the width of the curves. The 
main uncertainty in pinning down p c comes from the non- 
obvious and somewhat subjective extrapolation to t — ► oo. 



We should note that we calculated also P(t), the prob- 
ability that a cluster survives at least t steps (i.e. has 
"chemical radius" > t), the cluster size distribution P(n), 
and the spatial extent of clusters with n sites. All of them 
gave vastly more noisy signals (since we could not use a 
similar variance reduction trick as for M(t)) and were 
not used in estimating the critical point. 

Results for M(t) are shown in Figs.l to 3 for d = 6, 7, 
and 11. In all these figures we show results for bond per- 
colation. Results for site percolation are similar albeit 
somewhat more noisy. In the first two cases we checked 
explicitly that no cluster was larger than the virtual lat- 
tice size L (which was > 500 in both cases), so there are 
strictly no finite lattice size effects. For d > 9 this was 
no longer possible for the cluster sizes used here (typi- 
cally up to 10 4 — 10 6 sites), but we can easily convince 
ourselves that also there finite size effects are negligible. 

In each of the 3 figures the critical point p c is char- 
acterized by dM(t)/dt — ► for t — > oo. For d > 6 

we also have M(t) — > const for p = p c , while we see 
a logarithmic divergence in d = 6 as predicted by the 
renormalization group (l4| (see Fig.l). Unfortunately, 
the detailed behaviour of M(t) in d = 6 has not yet been 
calculated, though the results of 0] and the fact that 
v t (the exponent controlling the correlation time) is 1, 
suggest M(t) ~ [logi] 2 / 7 to leading order. Therefore, 
and since it is notoriously difficult to verify logarithmic 
terms (see, e.g., jl5[ O), we have not attempted any 
detailed analysis. 

From Figs. 2 and 3 we also see that corrections to scal- 
ing decrease strongly with dimension for d > 6. In Fig. 2 

we see a straight line for p — p c when plotting M(t) 
against 1/Vt, showing that the leading correction term is 
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in d — 7. Similarly, a straight line is obtained for 
d = 11 when using i~ 2 - 18 (Fig. 3). All these (and similar 
results for other values of d > 6 and for site percolation, 



p = 0.0786757 
p = 0.0786752 
p = 0.0786747 



1.23 
1.22 
1.21 
1.2 

^ 1.19 
1.18 
1.17 
1.16 
1.15 

0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 

1/t 1/2 

FIG. 2: Plot of M(t) versus t~ 1/2 for bond percolation in 
d — 7. Statistical errors are always smaller than half the 
distances between neighboring curves. 
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FIG. 3: Plot of M(t) versus t~ 



for bond percolation in 



d = 11. Statistical errors are again smaller than half the 
distances between neighboring curves. The exponent 2.18 is 
chosen since it gives the straightest line. 



not shown here) strongly suggest anomalous scaling 

M{t) = Moo - const/r (d) , (6) 

similar to the scaling for self avoiding walks in d > d c 
found in |17|] . But while the exponents were simply (d — 
d c )/2 in Wq, they seem to depend less trivially on d 
in the present case - although we cannot exclude the 
possibility that ui(d) = (d — d c )/2 also here, and the 
observed deviations are due to higher order corrections. 
The latter is indeed suggested by the results of |H . 

The constants M m defined in Eq.(^|) are plotted in 
Fig. 4 against d — 6 on doubly logarithmic scale. They 
seem to fall on parallel straight lines, suggesting a uni- 
versal law Moo - 1 ~ (d - 6)~ a with a = 0.73 ± 0.03. 
But a closer look reveals that deviations from this are 
significant (although they are small), suggesting that it 
holds neither for d — > oo nor for d — > 6 exactly. 
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FIG. 4: Log-log plot of M x — 1 against d — d c . Statistical 
errors are smaller than the data symbols. 
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FIG. 5: Log-log plot of the discrepancies between the simula- 
tion results and Eqs.(Q) and (^). Error bars are smaller than 
the sizes of the symbols. 



Let us finally discuss the p c values given in Table 1. 
They should be compared to the predictions 



Pcbond = a + 5s 3 /2 + 15s 4 /2 + 57s 5 



(7) 



and M 



Pcsite = s + 3s 2 /2 + 15s 3 /4 + 83.s 4 /4 + • ■ ■ (8) 



with s — l/(2d— 1). The dots in these equations stand for 
higher powers of s. It was suggested in [19| that they can 
be approximated, for site percolation at least, by adding 
2/3 of the last term, 

p c ,sitc ~ s + 3s 2 /2 + 15s 3 / 4 + 415s 4 /12. (9) 



The full series are presumably only asymptotic. It is thus 
a priori not clear whether any of these equations should 
be good approximations to the present data. From Ta- 
ble 1 we see that Eq.(^) is excellent in the range studied 
here, but it has wrong asymptotic behaviour and should 
be worse than Eq.(||) for d > 15. As seen from Fig. 5, 
the agreement with Eqs.(|7]) and (||) is indeed better than 
could have been expected: For bond percolation the dif- 
ference decreases roughly as s 71 (instead of s 6 ), while for 
site percolation it decreases as s 5 7 instead of s 5 . Obvi- 
ously the next terms in Eqs.(Q) and @ would be needed 
for a more detailed comparison. 

Finally we should remind of several heuristic formu- 
las for p c values on various lattices. All early ansatzes of 
this type were already refuted in ^TJ] because they contra- 
dicted Eqs.(0) or (||). More recently, such heuristics have 
been discussed again in pl| and in the papers quoted 
there. We have not attempted any detailed comparison 
in view of their complete lack of theoretical basis. 

In summary, we have presented vastly improved es- 
timates for percolation thresholds on high-dimensional 
hypercubic lattices. They should be compared to im- 
proved series expansions and/or rigorous bounds. At 
present such results are not available, partly because it 
had seemed that they could not be compared to any nu- 
merical estimates. Apart from this, the methods used 
in the present paper should also be of use in other sim- 
ilar problems. These include simulations of percolation 
backbones, conductivity exponents, percolation on more 
exotic lattices, directed percolation in high dimensions, 
and self avoiding walks. In all these cases both the hash- 
ing and the variance reduction should be of help in sim- 
ulating larger systems with higher precision. 

I thank Walter Nadler and Hsiao-Ping Hsu for discus- 
sions and for carefully reading the manuscript. 
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